# Run summary stats for survey analysis


# Setup ----
library(tidyverse)
library(xtable)

# Load survey data as "a"
# SOURCE: please directly contact the authors of You and Kousky (2024)




################################################################################
# Table 7: Summary Statistics, Survey Sample
################################################################################


# summary stats function -------------------------------------------------------
sum_stats = function(dataname, var1){                                           
  # Summary stats table to be passed to dplyr functions
  # dataname is a data frame
  # y is a column in that data frame
  var1 = enquo(var1)
  
  b0 = dataname %>% select(!! var1)
  
  b1 = dataname %>% filter(AnyIns==0 & !is.na(!!var1) & !is.infinite(!!var1)) %>% select(!!var1)
  b2 = dataname %>% filter(AnyIns==1 & !is.na(!!var1) & !is.infinite(!!var1)) %>% select(!!var1)
  
  a0 <- b0 %>% 
    summarise(Mean=format(round(mean(!!var1, na.rm = TRUE), 2), nsmall = 2), 
              Median=format(round(median(!!var1, na.rm = TRUE), 2), nsmall = 2),
              SD=format(round(sd(!!var1, na.rm = TRUE), 2), nsmall = 2)) %>% 
    data.frame()
  
  a0_i <- b1 %>% 
    summarise(NoIns=format(round(mean(!!var1, na.rm = TRUE), 2), nsmall = 2)) %>% 
    data.frame()
  
  a0_ii <- b2 %>% 
    summarise(AnyIns=format(round(mean(!!var1, na.rm = TRUE), 2), nsmall = 2)) %>% 
    data.frame()
  
  t1 <- t.test(b1, b2)
  t_stat <- paste0(format(round(t1$statistic, 2), nsmall = 2), 
               ifelse(t1$p.value<=0.01,"***",
                      ifelse(t1$p.value<=0.05,"**",
                             ifelse(t1$p.value<=0.1,"*",""))))
  
  a3 <- cbind.data.frame(a0, "", a0_i, a0_ii, t_stat)
  rownames(a3) <- colnames(b0)[1]

  
  return(a3)
}



# disaster experience -----------------------------------
a = a %>% data.frame()

m1  = sum_stats(a, enough_money_low)
m2  = sum_stats(a, financial_yr_post_veryworse_ind)
m3  = sum_stats(a, Unmet_Needs)
m4  = sum_stats(a, damage_repair_time_cat)

m5a = sum_stats(a, damage_home_Yes)
m5b = sum_stats(a, damage_contents_Yes)
m5c = sum_stats(a, service_disrupt_Yes)
m5d = sum_stats(a, evac_decision_Yes)

m6a = sum_stats(a, damage_severity_home_num)
m6b = sum_stats(a, damage_severity_contents_num)
m6c = sum_stats(a, service_disrupt_cost_num)
m6d = sum_stats(a, evac_totaldollars_000)
m6e = sum_stats(a, damage_car_Yes)
m6f = sum_stats(a, lostincome_dummy)
m6g = sum_stats(a, costs_additional_Yes)






# home characteristics ----------------------------

n1  = sum_stats(a, home_residents_num_new)
n2  = sum_stats(a, children_disability_elderly_pets)
n3  = sum_stats(a, home_tenure)
n4  = sum_stats(a, home_mortgage_dummy)
n5  = sum_stats(a, renter)
n6  = sum_stats(a, single_family)





# household characteristics ----------------------------
n0  = sum_stats(a, savings_yes)
n5a = sum_stats(a, hhi_low)
n5b = sum_stats(a, hhi_median)
n5c = sum_stats(a, hhi_high)

n6a = sum_stats(a, employment_fulltime)
n6b = sum_stats(a, employment_parttime_selfemployed)
n6c = sum_stats(a, employment_retired)
n6d = sum_stats(a, employment_Other)

n7  = sum_stats(a, nonwhite)






summary_table = rbind(
                      m1,  m2,  m3,  m4,  
                      m5a, m5b, m5c, m5d,
                      m6a, m6b, m6c, m6d, m6e, m6f, m6g,
                      n1,  n2,  n3,  n4,  n5, n6,
                      n0,  n5a, n5b, n5c, 
                      n6a, n6b, n6c, n6d, n7
                      )



xtable(summary_table)





